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An experimental, computational and theoretical investigation was carried out to study 
the aerodynamic loads acting on a relatively thick NACA 0018 airfoil when subjected to 
pitching and surging, individually and synchronously. Both pre-stall and post-stall angles of 
attack were considered. Experiments were carried out in a dedicated unsteady wind tunnel, 
with large surge amplitudes, and airfoil loads were estimated by means of unsteady surface 
mounted pressure measurements. Theoretical predictions were based on Theodorsen’s and 
Isaacs’ results as well as on the relatively recent generalizations of van der Wall. Both 
two- and three-dimensional computations were performed on structured grids employing 
unsteady Reynolds-averaged Navier-Stokes (URANS). For pure surging at pre-stall an- 
gles of attack, the correspondence between experiments and theory was satisfactory; this 
served as a validation of Isaacs theory. Discrepancies were traced to dynamic trailing-edge 
separation, even at low angles of attack. Excellent correspondence was found between 
experiments and theory for airfoil pitching as well as combined pitching and surging; the 
latter appears to be the first clear validation of van der Wall’s theoretical results. Although 
qualitatively similar to experiment at low angles of attack, two-dimensional URANS com- 
putations yielded notable errors in the unsteady load effects of pitching, surging and their 
synchronous combination. The main reason is believed to be that the URANS equations do 
not resolve wake vorticity (explicitly modeled in the theory) or the resulting rolled-up un- 
steady flow structures because high values of eddy viscosity tend to “smear” the wake. At 
post-stall angles, three-dimensional computations illustrated the importance of modeling 
the tunnel side walls. 
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Nomenclature 


a 

Angle of attack, ° 

Ad 

Angle of attack oscillation amplitude, ° 

7 

Vortex sheet strength, m/s 

A 

Wave length of the oscillating, m 

P 

Air density, kg/m 3 

a 

Relative velocity amplitude of the oscillating 
free stream 

<t> 

Phase angle, ° 

LU 

Angular frequency, 1/s 

OJ s 

Rotational speed, rad/s 

r 

Circulation, m 2 /s 

a 

Non-dimensional pitch axis 

c 

Airfoil chord length, m 

k 

Reduced frequency, uic/ (2 u s ) 

l 

Coefficients from Isaacs 

m 

Wave number 

p 

Static pressure, Pa 

t 

Airfoil thickness, m 

u 

Free stream velocity, m/s 

Ci 

Lift coefficient 

c p 

Pressure coefficient 

C(k ) 

Theodorsen function 

m 

Real part of the Theodorsen function 

G(k) 

Imaginary part of the Theodorsen function 

H 

Coefficient from van der Wall 

J 

Bessel function of the first kind 

L 

Lift, N 

R 

Radius of a vertical axis wind turbine, m 

Re 

Reynolds number 

Subscript 



qs 

Quasi steady 

ref 

Reference 

s 

Steady 

st 

static 

us 

Per unit span 



Operators 



0 

Time averaged 

0 

Phase averaged signal 

/ 

Turbulence of a signal 

A 

Difference 


I. Introduction 

During the first half of the 20th century, the study of unsteady aerodynamics was motivated by problems 
associated with wing flutter, the estimation of helicopter blade loads, and the effect of wind gusts on aircraft. 
Many of these problems remain relevant today. In recent years, the stability and survivability of small low 
flying air vehicles that are exposed to highly unsteady winds have produced renewed interest in the topic. 
Furthermore, the desire for robust design of wind turbines - whose blades are exposed to highly unsteady 
flows produced by, inter alia, yaw misalignments and atmospheric turbulence - requires an urgent need for 
the accurate prediction of unsteady loads. 1 In this paper, two sources of unsteadiness are explored. On the 
one hand, the attached flow in unsteady conditions like an oscillating free stream or a pitching airfoil at 
low angles of attack is studied. Although the flow remains fully attached, significant unsteady lift overshoot 
has been documented, see Refs. 2-4. On the other hand, the phenomenon of dynamic stall is investigated, 
which is characterized by the shedding of a strong vortex across the suction surface that induces rapid load 
excursions. Dynamic stall itself is already a highly unsteady phenomenon; nevertheless, combined studies 
including an oscillating free stream are rare and a lack of validation still exists. 

Considering the fully attached airfoil, theoretical approaches to unsteady aerodynamics commenced about 
two decades after the work of Kutta 5 and Joukowski (Zuhkouski). 6 Calculations of forces acting on a flat 
plate in unsteady plunge motion are published in Wagner. 7 It also provided a first approach to calculate the 
lift of airfoils in pitching motion, which is the basis for the development of a complete theory for pitching 
motion in Glauert. 8 Challenged by the problem of wing flutter, a landmark innovation occurred with the 
development of a general analytical solution for airfoils undergoing angle of attack oscillations and plunge 
motion in Theodorsen. 2 All velocity potentials were calculated for a steady stream potential flow and the 
unsteady circulation was computed, where the wake vorticity is determined by the Kutta condition. Today, 
we refer to the Theodorsen function that includes multiples of the reduced frequency as its argument. 

The need for more accurate estimations of helicopter blade loads motivated Isaacs 3 to extend Theodorsen’s 
work to include an oscillating free stream with no phase lag between the leading and trailing edges. A year 
later, the theory was extended to include angle of attack oscillations at the same frequency. 9 Greenberg 10 
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developed a solution for the dynamic lift and moment of a flat plate in an oscillating free stream undergoing 
angle of attack oscillations and plunge motion. Furthermore, this approach included an arbitrary phase 
lag between the oscillating free stream and the periodically pitching wing. However, this modulation of 
the unsteady free stream also implied a fore-aft wing motion. More than 30 years later, a new theory was 
developed, which was directly motivated by helicopter aerodynamics. Its focus was the prediction of the 
dynamic drag behavior under unsteady free stream. The derivation was based on an order of magnitude 
approximation to the in-plane perturbation velocity. 11 Its author asserted that this assumption is more 
plausible than the high frequency assumption in Ref. 10. However, this theory was limited to small free stream 
fluctuations. In 1991, an extensive review of existing theoretical approaches was provided and solutions for 
harmonic plunge motion and unsteady angle of attack variations including arbitrary multiples of the free 
stream harmonic were provided, see Refs. 4, 12-14. The review concluded that Isaacs’ theory is the only 
“exact theory” without additional simplifications. Furthermore, significant deviations of Greenberg’s theory 
for relative velocity amplitudes higher than a = 0.4 were determined. Based on Theodorsen’s approach, an 
analytical framework was developed for variable geometry airfoils including camberline elasticity in arbitrary 
motion. However, this ansatz did not consider the effect of an unsteady free stream, see Ref. 15. 

In addition to fully attached flow configurations, detached flow phenomena are investigated in this study 
as well. Dynamic stall is a fascinating fluid dynamics phenomenon observed on rotorcraft and wind turbines 
when their blades pitch rapidly beyond the static stall angle, a s t . Excellent reviews of dynamic stall and its 
prediction have been published. 16,17 Carr 16 describes approximately eleven stages of airfoil dynamic stall at 
high Reynolds numbers Re > 10 6 . In summary: reversed flow initially appears on the airfoil upper surface 
and spreads; the dynamic stall vortex (DSV) forms near the leading edge producing increases in lift and 
moment; finally, the vortex sheds from the airfoil leading to lift stall. Ekaterinaris et al. 17 also provided a 
summary of many CFD efforts to predict dynamic stall, and describe important issues related to turbulence 
and transition modeling when using unsteady Reynolds averaged Navier-Stokes (URANS) methodologies. 
Many other CFD studies for dynamic stall have been conducted since the time of the review, including flow 
control, low Reynolds number studies, and three-dimensional applications (e.g., Refs. 18-23). 

On wind turbines, dynamic stall manifests in different ways depending upon the type of turbine. On 
vertical axis wind turbines (VAWTs), blades stall dynamically at low blade tip speed to wind speed ratios, 
A = oj s R/u s (oj s is the rotational speed, R is the turbine radius and u s is the wind speed), resulting in loss 
of power and unsteady loads. Understanding and modeling of dynamic stall, particularly due to the current 
renewed interest in VAWT power, is of major scientific and technological importance (e.g., Refs. 24,25 and 
many others). For structural reasons, these blades are usually thick ( t/c ~ 18%) and the Reynolds numbers 
can be low, typically between 5 • 10 4 < Re < 10 6 , particularly on smaller urban and off-grid machines. 

Experiments involving an unsteady free stream have been achieved using many novel designs (see Refs. 26- 
30 and the references therein). In experiments involving stalled flows, some remarkable results have been 
obtained when the free stream varies according to u(t) = it s (l + a cos (u/t)), where er is the relative velocity 
oscillation amplitude, and the angular frequency, ui = 2irf. For example, Gursul and Ho 31 measured mean 
lift coefficients on a NACA 0012 airfoil at a = 20° far in excess of comparable quasi static values. The lift 
depended strongly on velocity oscillation amplitude (also see Ref. 32) and frequency, producing instantaneous 
lift coefficient values that were, amazingly, larger than 10, with optimum reduced frequencies k = luc/2ii s « 
0(1). The source of the high lift is a large vortex that is momentarily trapped on the airfoil and then swept 
downstream. Similar observations were made on delta wings at high angles of attack (Refs. 33-36). The 
only investigation of two-dimensional dynamic stall combined with an unsteady free stream appears to be 
that of Pierce et al. 37 Although their angle of attack range was small a(t) = ±4°, the authors observed a 
large variation in the post-stall pitching moment coefficient as a result of the flow unsteadiness. 

As discussed above, fully attached flow in unsteady conditions yield significant lift overshoots. Thus, it 
is surprising that many of these theories, in particular those of Greenberg and Isaacs, have not been fully 
validated experimentally or numerically, especially for large free stream oscillation amplitudes. 38 Favier et 
al. 39,40 investigated a pitching NACA 0012 airfoil in unsteady free stream. The wind tunnel generated high 
reduced frequencies up to k = 1.6 and relative velocity amplitudes of more than er = 0.35. Although the 
airfoil lift shows significant dynamic effects, a comparison to Isaacs’ theory was not considered. More recently, 
Granlund et al. 41 investigated a NACA 0009 experimentally over a broad range of reduced frequencies at 
small relative velocity amplitude of a = 0.1. The reason for this lack of validation is unclear, but apparently 
the existing experimental facilities lack the large amplitude unsteady parameter range. Various wind tunnel 
challenges include fan stall, and large lag effects due to inertia and acoustic resonance. Indeed, a review 
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of facilities reveals that the vast majority of them simulate only changes in airfoil angle of attack, namely 
pitching or plunging. 42-44 Tunnels that produce an unsteady free stream are rare. The most common 
approach is to modify a standard steady wind tunnel to produce unsteady flows, 37 ’ 45-47 typically using 
rotating vanes, louvers or some variation of this idea. Some tunnels combine the independent capabilities 
of angle of attack and wind speed variation. 39, 48,49 An exception is the approach of Favier et al., 40 which 
employs a two degrees of freedom construction to pitch and translate the airfoil. Recently, an unsteady 
wind tunnel was developed to produce large amplitude oscillations of the free stream in Furman et al. 50 
Problems of fan stall, large inertial effects, and acoustic resonance were overcome during the initial design 
and testing phases. The tunnel proved to be ideally suited to validate large amplitude unsteady effects and, 
in particular, to assess the validity of theoretical approaches. 

On helicopter blades, on VAWTs, and on HAWTs exposed to yaw misalignment, dynamic stall is ac- 
companied by other sources of unsteadiness. In particular, the physical reality involves a dynamic change 
in pitch angle with a simultaneous change in flow speed. However, this simultaneous variation/oscillation 
is almost never reproduced experimentally in the context of airfoil dynamic stall. In the field of VAWTs, 
a notable exception is the study by Ferreira et al., 24 who made two-dimensional particle image velocime- 
try (PIV) measurements on a small VAWT at Re < 10 5 and tracked the dynamic stall vortex at different 
phase angles, <f> = tot. However, we have chosen a NACA 0018 airfoil as the focus of our study. 51 Most 
of the static testing of NACA 0018 airfoils has been conducted at high Reynolds numbers appropriate to 
the aircraft industry and large turbines. 52,53 Less research has been conducted at Re < 10 6 . Aerodynamic 
performance is significantly affected as the Reynolds number decreases below Re = 3 • 10 5 as shown in the 
investigations of Jacobs and Sherman 54 and Raghunathan and Ombaka, 55 and more recently in Timmer 56 
and Gerakopulos et al. 5 ' Similarly, dynamic stall on NACA 0018 airfoils has been studied mainly at rela- 
tively large Reynolds numbers, Re > 10 6 (see Refs. 58-60, where the emphasis is on modeling stall and flow 
reattachment. Dynamic stall studies on NACA 0018 airfoils at Re < 3 • 10 5 are unusual 55 ). 

The first objective of this investigation is to describe a combined theoretical, numerical and experimental 
study aimed at validating the theories of Greenberg and Isaacs for large amplitude free stream oscillations 
of cr = 0.5. Three different theoretical approaches to predict the overshoot are compared and the ability 
of URANS to capture the unsteady flow physics is assessed. For this purpose, the lift, the leading-edge 
pitching moment, and the pressure distribution along the chord are valuable parameters for comparison. The 
second objective is to experimentally and computationally study dynamic stall under idealized conditions 
representative of wind turbine rotor blades. Dynamic stall occurs suddenly and is influenced by the pressure 
gradient along the chord, the boundary layer, the global dynamics in the system (expressed by the reduced 
frequency), etc. This is a challenging setup for CFD and the conducted experiments represent a great basis 
for validation. 

II. Description of CFD code used for unsteady flow predictions 

Alongside the experimental effort and Isaacs’ theory, a companion CFD study is being conducted for 
the same configurations. At this time, two-dimensional steady RANS and unsteady URANS calculations 
have been employed for several different cases. The latter computations include simulation of the unsteady 
free stream and moving airfoil, and are compared with the experimental data and theory. Also, a three- 
dimensional steady RANS calculation including the test section (side walls, floor and ceiling) at several 
angles of attack has been performed in order to explore the three-dimensionality of the flow at higher lift 
conditions. Due to the higher computational cost of the three-dimensional calculations, this has not yet been 
extended to unsteady flows. 

The structured grid CFD code CFL3D 61 is employed for the RANS and URANS computations. It is 
an upwind finite volume code that solves the full Navier-Stokes equations. CFL3D can solve multiple zone 
structured grids that are connected in an one-to-one, patched, or overset manner. It possesses the Spalart- 
Allmaras (SA), 62 Menter SST, 63 and SSG/LRR-RSM-2012 second moment Reynolds stress 64 turbulence 
models (among others) to simulate turbulent flows. For all results to be shown in this paper, only the SA 
model was employed. It was run as fully turbulent, with no presumed transition on the airfoil. CFL3D 
is globally second order spatially accurate, and can be used to solve time accurate problems with moving 
bodies and with flow control applied at walls. 65 

For time accurate problems, CFL3D utilizes pseudo time stepping and achieves second order temporal 
accuracy. With pseudo time stepping, sub-iterations are used to reduce the linearization and factorization 
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errors, and advance the solution in pseudo time to the next physical time. 

For time accurate URANS computations, the decomposition for forced periodic cases can be written as 

F(x,t) = (F(x,t)) +F'(x,t) (1) 

where F(x,t) is the total signal, the ( F(x,t )) term is the phase average, and F'(x,t ) is the turbulence. This 
form resembles standard Reynolds decomposition, except that the flow is split into a slowly varying mean 
(phase average) and a random fluctuating part. As a result, the final conservation equations in terms of 
phase averaged variables ( F(x,t )) are identical in form to the standard Reynolds averaged equations. In 
other words, when URANS is used for a forced periodic flow field such as that produced by oscillations of 
the airfoil or free stream perturbations, the computation (depending on the turbulence model and the case) 
often eventually settles down to a nearly exactly repeatable periodic variation. If and when this repeatability 
is attained, any single point during the cycle corresponds to a phase averaged result from the experiment. 



Figure 1. Part of the tunnel grid with the NACA 0018 airfoil used by CFD (zone boundaries are shown in 
bold red lines). 


In order to include the varying free stream speed in the CFD, at the outflow of the tunnel grid the 
boundary condition sets a sinusoidal time varying back pressure, resulting in variable speed through the 
tunnel. This strategy attempts to closely mimic the tunnel conditions. However, several difficulties have 
been noted with this approach. First, it is a trial-and-error process to discover the back pressure mean 
magnitude p and change A p that will yield a desired tunnel velocity variation. Second, the conditions 
achieved are a strong function not only of the back pressure, but also of the frequency of oscillation. In 
other words, at one frequency a certain p and A p will produce the desired tunnel velocity variation, but at a 
different frequency the same numbers no longer work, and a new set must be found by trial-and-error. And 
finally, when performing simultaneous tunnel speed variation and airfoil pitching, the phase of the tunnel 
back pressure must be adjusted so that the tunnel velocity is phased appropriately with the airfoil pitching. 

For simulating the airfoil pitching, circular grid zones were constructed around the airfoil, and these were 
dynamically rotated during the computation. In this process, information is passed to/from the surrounding 
tunnel grid via patched interfaces. A picture of the fine grid is shown in figure 1, focusing on the region near 
and just upstream of the airfoil. The two zones with circular outer boundary are rotated to a = 2° in this 
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view. The grid includes the tunnel’s upstream settling chamber and contraction; the inflow in the grid is at 
x = —4.5m, and the contraction ratio is 7.91. The outflow location is at x = 3.77m. Note in the CFD that 
z is “up” and y is spanwise. 

The boundary conditions at the airfoil surface are viscous (no-slip) adiabatic. The tunnel top and bottom 
walls are simulated via inviscicl slip-wall conditions. At the tunnel outflow, the pressure ( p/p re f ) is specified 
(either constant or sinusoidally varying, depending on the case), while all other quantities are extrapolated 
from the interior of the domain. At the tunnel inflow, the total pressure and total temperature are specified 
at free stream (reference) levels, and the Riemann invariant is extrapolated from the interior of the domain. 


III. Experimental replication of blade conditions 


The wind tunnel is driven by a 75kW blower, has an 8:1 contraction ratio and a top wind speed of 
55 m/s. 50,51 The wind tunnel was designed specifically to accommodate smooth flow speed variations by 
selecting a backward bladed impeller that operates smoothly even when fully stalled. Flow distortion in the 
upstream test section and wind tunnel turbulence are less than 0.2% and 0.1% respectively. The test section 
is equipped with large rings that incorporate transparent Plexiglas® windows (see figure 2). The rings are 
driven by a servomotor via a series of two steel reinforced belt drives. In addition, the floor and ceiling are 
also made of Plexiglas® rendering full optical access to the test section. In this study, the NACA 0018 airfoil 
is attached to, and mounted between, the Plexiglas® walls. The rings are able to perform arbitrary pitching 
motion of the airfoil through a full 360° at maximum rate of 150°/s. This custom built test section with 
arbitrarily large angle of attack range, combined with high pitching speeds and full optical access renders this 
setup unique. The system was programmed using Lab VIEW™ with the so-called CompactRIO™, which 
is a realtime industrial controller made by National Instruments. Considerable expense and manpower has 
been invested in perfecting and validating this test section. Design of the system took into account the 
stiffness of the belts and motor response to the loads. At peak loads, the maximum error was calculated to 
be < 0.25°; direct measurements confirmed the fidelity of the design with errors never exceeding 0.2°. 


servo motor 



Figure 2. View of the test section showing the servomotor and belt drives as well as the approximate airfoil 
location. 

The well-known method of employing louver vanes was adopted for the present tunnel. The vanes 
dynamically vary the cross sectional exit area of the tunnel and thus, vary the head-loss resulting in unsteady 
velocity oscillations. The present implementation in figure 3 involves a louver mechanism with 13 counter 
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rotating vanes driven by a 0.75kW servo motor with a 5:1 gear ratio, coupled with PLC controller. This 
enables the dynamic control of the vanes’ position, and the programming of various arbitrary rotation 
profiles, thus creating different temporal velocity profiles inside the upstream test section. A comparison of 
the measured, sinusoidal velocity profile with an ideal sine yield a maximum deviation of 2% at a velocity ratio 
of 0.5. 66 Thus, the entire free stream velocity profile is expected to be sufficient for dynamic measurements. 



gear 1:5 


servo motor 


cog wheels 


vanes 



Figure 3. Photograph of the 13 louver vans at the end of the test section in open and closed position. 


The lift is calculated by means of 40 pressure taps located on the centerplane of the model. The measured 
static pressure, which acts normal to the surface, is weighted by the half distance to the neighboring pressure 
taps and transformed in the coordinate system of the wing chord. The summation yields the lift and the 
pressure drag (drag due to friction is not quantified). The cross product of the static pressure at each pressure 
tap and the distance to the leading edge gives the leading edge pitching moment. The phase reconstruction 
is based on the averaged free stream velocity of the two hotwires. Taking into account that the amplitude 
of the free stream oscillation varies slightly, each single period was fitted by an ideal sine to avoid any 
unphysical scatter in the data. Each measurement consists of at least 150 periods. The data are averaged 
at each a s t ep = 0.5° with a window size of ±0.3° and at each <f> s te P = 2° with a window size of ±1°. 

IV. Theoretical unsteady lift prediction in pre-stall configuration 

As discussed above, unsteady conditions may yield significant unsteady lift effects on airfoils even when 
the flow remains fully attached. For this scenario, several theories predict the lift response assuming a 
flat plate in incompressible potential flow. 2-4,9,10,13 - 14 The assumed inflow and the airfoil geometry are 
homogenous and two-dimensional. The infinite wake itself is planar, starts at the trailing edge and no 
diffusion is considered. The wake convects downstream along the x-axis without any vertical or horizontal 
displacement due to self induction. The amplitude of the oscillating free stream is limited to values er < 1 to 
avoid reverse flow. The assumption of potential flow requires fully attached flow in experiments for a valid 
comparison. Thus, trailing edge stall or laminar separation bubbles lead to significant deviations. 

For a pure harmonic pitching motion in constant free stream, Theodorsen 2 developed a closed form 
solution to predict the unsteady lift overshoot. It includes the so called Theodorsen functions C(nk) = 
F(nk) + iG(nk ), which represent the basis for all further solutions. The unsteady lift generated by a 
sinusoidal angle of attack oscillation a{(j>) = ao + /3 q sin(<5) is given in equation (2). The steady lift is defined 
as Lq = Trpcu^ao- 

= 27r(0.5fc(/? o cos(</>) + ak/3o sin((/>)) + ao + Po sin (cj))F{nk) 
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fc(0.5 — a)f3o sin (<j))G(nk) + /3 q cos ((j))G(nk) + k( 0.5 — a)/ 3 q cos ((j>)F(nk)) 


(2) 


Considering an airfoil at a constant angle of attack exposed to an unsteady free stream, Greenberg 10 
and Isaacs 3 provided solutions to predict the unsteady lift response. However, for large relative velocity 
amplitudes of a > 0.4, Greenberg’s approach yields significant deviations from theory 4 and experiments. 66 
Thus, only Isaacs’ results are discussed here briefly and shown in equation (3), assuming a sinusoidal free 
stream u(<j>) = tt s (l + a sin(^)). The lift coefficient from this equation is used instead of the lift coefficient 
from ( 1 -|_ <TS in( u ,t ))2 ■ This interpretation is independent of the dynamic pressure variations of the 

free stream and hence represents the pure unsteady lift response. 


Ci(t) 

Cl,qs 


— . . [1 + 0.5a 2 + a(l + 9f(Zi) + 0.5a 2 ) sin (wt) 

(1 + asm(uA)U 

OO 

+a(5i(Zi) + 0.5fc) cos (cut) + a ^ (SR(Z m ) cos (muit) + 3(Z m ) sin(mwf))] 

m—2 


(3) 


Note that Z m includes an infinite summation, and J represents Bessel functions of the first kind. 


OO 

lm = TYl[ Z) ^ ^ [Tn(d n -)_ m ( Ua) T n _ m (?7.a)) ~t“ iCjn ( J n ^rn (ua ) + T n _ m (ua))] (4) 

n — 1 


F n 1 _ J n+ i(na) - J n -i(na) j F(nk) 

G n j n 2 | G(nk) 

Isaacs was the first to publish a closed form solution for the combined case of a sinusoidal free stream 
oscillation and a harmonically pitching wing. 9 Based on this approach, van der Wall generalized Isaacs’ 
ansatz and included plunge motion and an arbitrary location of the pitch axis. 4 The unsteady lift response 
for a sinusoidal free stream u(<j)) = u s (l+asin(<^>)) and an in-phase pitching motion a(cj>) = ao(do+ai s sin(</>)) 
is given in equation (6). 


cm 

Cl,qs 


OO 

(1 + asin(</>))~ 2 [((l + 0.5 a 2 )a 0 + aa i s )(l + asin(^)) + ^ (Z m cos(m</>) + l' m sin(m0)) 

m—1 


+0.5/c((aa o + &is) cos(^) + kaa i s sin(^) + aai s sin(</>))] 


( 6 ) 


Due to the combined occurrence of the two degrees of freedom, the definition of F n + iG n changes and 
includes more parameters. 4 


F n + iG n 

H n 

K 


n 2 [F(nk) + iG(nk))(H n + iH ' n ) 

0.5(J„ + i - J n _i)(aa 0 - ai a ) - (na) _1 (2 J n a ls ) 
o’ -1 J„(-fc(0.5 - o)ai a ) 


(7) 

(8) 
(9) 


V. Results 

In a first step, the baseline of the NACA 0018 experiments and the two- and three-dimensional CFD 
simulations are described. Then, the unsteady lift response in fully attached flow at low angles of attack 
are presented. An oscillating free stream, pitching motion and the combination of both unsteady effects are 
compared to theoretical predictions. Finally, implications for the airfoil pitched beyond the static stall angle 
are discussed. 

In all CFD cases, both medium and fine grids were tested, and the differences in results were found to 
be minimal. Unless otherwise noted, fine grid results are shown for all cases. For unsteady flows, the CFD 
also explored the effects of number of subiterations (per physical time step) and time step value itself on the 
results. Generally, use of at least 10 subiterations in combination with a time step yielding at least 600 steps 
per period were found to be sufficient. Unsteady computations were initiated from steady-state solutions, 
and continued until a repeatable periodic variation was reached (typically after 2-3 cycles). 
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A. Steady results from CFD and experiment 


In order to assess the angle of attack range over which CFD is likely to be reliable, steady CFD results 
were obtained over a range of angles of attack and compared with experimental data. Results are shown in 
figure 4. Two different grid sizes were employed to assess the influence of discretization errors. The fine grid 
(f) had 97,280 grid cells, and the medium grid (m) had 24,320 grid cells. As seen in figure 4, there was little 
influence of grid size on the resulting lift coefficient. 



Figure 4. Steady lift coefficients of a NACA 0018 elaborated by quasi-steady experiments (continuous pitch 
sweeps) and two-dimensional steady CFD at Re = 3 ■ 10 5 . 

In the experiment, angles of attack well beyond stall were considered, and hysteresis was noted in the 
range 15° < a < 21°, as seen in figure 4. Data were acquired using continuous pitch sweeps at very low 
k = 0.000085. The maximum lift coefficient Cl (near a = 14.5°) was around 1.1. In CFD, results started 
to deviate from the experiment above a of 8-10°. Near here, the lift coefficient slope of the experiments 
was reduced significantly below the theoretical slope (assuming potential flow) of 27 t, suggesting the onset of 
trailing edge stall. At a = 14.5°, the computed Cl of about 1.3 significantly overpredicted the experimental 
value. There is no doubt that two-dimensional CFD is incapable of matching the experiment at higher angles 
of attack, so the two-dimensional steady computations were taken no further than 14.5°. 

Two three-dimensional steady RANS computation were run at a = 2° and 14.5°, with viscous tunnel side 
walls included. For these computations, the two-dimensional fine grid was extended in the spanwise direction 
using 64 grid cells over y = 0.305m (half the full tunnel width), for a total grid size of 6.2 million grid cells. 
The grid was clustered at the side wall, where viscous adiabatic boundary conditions were imposed. At the 
centerplane, symmetry boundary conditions were imposed. 

Surface pressure coefficient contours along with surface streamline patterns are shown in figures 5 and 6. 
At a = 2°, the flow is attached and two-dimensional over most of the wing, with only a very small region of 
three-dimensionality very near the side wall. At the centerplane, the flow was fully attached, satisfying the 
Kutta condition. This enables a comparison to the previously discussed theories that assume a flat plate in 
potential flow. Hence, this setup at low angles of attack provides a unique possibility of validating CFD, 
experiments and theory. Integration of the CFD surface pressure coefficients along the centerplane of the 
3-D computation yield a lift coefficient value per unit span of approximately 0.203, which is very close to 
experiment and to the 2-D computational result shown in figure 4. 

At the higher angle of attack of a = 14.5°, there is clearly a substantial influence of the side wall, as 
seen in figure 6. Note that this computation predicts significant regions of separated flow, and is not fully 
converged; therefore, this result should be viewed only qualitatively. Even though the flow appears to be 
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Figure 5. Surface pressure coefficient contours and surface streamlines for three-dimensional steady compu- 
tation on fine grid at a = 2°. 


reasonably two-dimensional over a large portion of the central region of the wing, the three-dimensionality 
has still affected the flow in that region. This can best be seen in a plot of surface pressure coefficient taken 
at the center plane position. In figure 7, the three-dimensional computed C v is compared to results from 
the two-dimensional computation. The three-dimensional results show significantly lower suction values, in 
better agreement with experiment at the center plane location. Integration of the CFD surface pressure 
coefficients along the centerplane of the 3-D computation yield a lift coefficient value per unit span of 
approximately 1.17. This value is much closer to experiment than the 2-D computational result. 

A similar influence of three-dimensional sidewall structures on the assumed two-dimensional flow near 
the tunnel centerplane was noticed previously in Rumsey et al. 67 for a high lift multi-element airfoil. In other 
words, a local region of apparent two-dimensionality over a region of a wing in a tunnel is not a guarantee of 
two-dimensional flow in the purest sense of the word. Three-dimensional structures near the side walls still 
influence the entire flow field (typically the higher the angle of attack, the greater the influence). 

The experimental data of the quasi steady pitch up and pitch down do not show any appreciable deviations 
at a = 14.5° (figure 7); i.e., no hysteresis effects exist at this angle of attack. Furthermore, the measured 
pressure distribution suggests the possible presence of a laminar separation bubble at approximately x/c = 
0.07. Alternatively, the slight deviation in the surface pressure distribution may be attributed to a spanwise 
blowing slot (height 1.2mm) that can be used for active flow control. Even though the entire slot was sealed 
with thin adhesive tape, 51 the surface discontinuity may have altered the transition process. Furthermore, 
figures 6 and 7 both suggest trailing edge separation over a broad region of the airfoil. At the trailing edge 
itself, reverse flow occurs and the Kutta condition is violated. This renders any comparison to potential flow 
theory impossible at these high angles of attack. 
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Figure 6. Surface pressure coefficient contours and surface streamlines for three-dimensional steady compu- 
tation on fine grid at a = 14.5°. 



Figure 7. Steady surface pressure coefficients, NACA 0018, a = 14.5°. 
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B. Comparison of the unsteady lift overshoot of a NACA 0018 facing an oscillating free 
stream 

Validation of both the experimental and computational approaches in relatively simple fully attached flows 
was considered important. The first step was an investigation of the influence of an unsteady freestream on 
an airfoil at a fixed low angle of attack. A comparison to experiments and potential flow theory assesses the 
reliability of 2-D URANS. 

The results are presented in the form of lift coefficient ratio, as shown in figure 8. This kind of presentation 
compensates for the strong dynamic pressure variations during the free stream velocity variations; and hence, 
only the net unsteady effects are highlighted. Quasi steady experiments were performed by slowly varying the 
wind tunnel speed at k ~ 0( 10” 4 ). In the CFD, reproducing quasi steady conditions was computationally 
difficult because the boundary conditions needed to be determined iteratively at each time step, as described 
earlier. Hence, steady computations were performed instead at 22 corresponding free stream conditions, and 
then a 6th-order best-fit polynomial was used to create an estimate for the quasi steady CFD lift variation. 
The unsteady CFD runs, yielding C; as a function of <j>, were conducted at the appropriate k by time- varying 
the tunnel back pressure in order to achieve the desired tunnel velocity variation. 



<H°] 


Figure 8. Experiments, CFD, and theory for the unsteady to steady lift ratio as a function of phase angle <j> 
for a NACA 0018 at Re = 2.5 • 10 5 and cx. — 2° subjected to an unsteady free stream oscillation of a = 0.5 and 
k = 0.074. 

In this particular case, the NACA 0018 is exposed to a sinusoidal free stream oscillation with an amplitude 
of 50%, resulting in u((f>) = w s (l + 0.5 sin(</>)). The time averaged free velocity is u s = 11.1 m/s which 
corresponds to Re = 2.5 • 10 5 . The angle of attack was kept constant at a = 2° and the reduced frequency 
was k = 0.074. Over the first half of the cycle (when freestream velocity is high), the lift experiences a deficit 
compared to quasi steady values. Over the second half of the cycle (when freestream velocity is low), the 
trend reverses and the lift exhibits overshoots compared to quasi steady values. 

In figure 8, near the maximum lift overshoot, the experiments (shown as black dots) exhibited an oscil- 
lation at 260° < (f> < 320°. This phenomenon possessed a high frequency and the corresponding reduced 
frequency is approximately k = 0.7. By means of measured, phase averaged pressure distributions in this 
range, this was traced to the formation and shedding of a recirculation bubble directly upstream of the trail- 
ing edge. This high frequency phenomena affected the flow field at the trailing edge by modifying the Kutta 
condition. Thus, the circulation varies and the shed wake vortex sheet is modified as well to keep the global 
circulation constant in the system. Finally, the modified wake vortex distribution leads to another measured 
lift overshoot compared to the theory. It is not clear whether this is caused by low Reynolds number effects, 
the relatively thick airfoil, or their combination. The solid black line represents an approximation of the 
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experimental results by a Fourier series including only the first two harmonics. 

Regarding the CFD, although the trends are correct, the overall, or integrated, lift enhancing effect is 
underpredicted. To explain this, the results of the potential flow theories 9, 10 are also presented here. The 
overall lift enhancement is greater for the theory of Isaacs (blue line), given by equation (3), yielding a 
better agreement with the experiments. Nevertheless, both theories account for the unsteady shed vorticity 
at the trailing edge in the form of a semi-infinite wake. Greenberg made a high frequency simplification 
that is physically equivalent to assuming a quasi steady convection velocity of the shed wake. 14 Thus, for 
the relatively high amplitude oscillations of o = 0.5 considered here, Greenberg’s theory is considered less 
physically representative than Isaacs’. In any case, the modeling of the wake clearly plays a crucial role in 
predicting the unsteady lift. URANS, however, does not resolve the unsteady flow vorticity, or structures 
resulting from rollup of the vortices in the wake, because the high values of eddy viscosity tend to “smear” 
all wake structures. 

For example, figure 9 is a snapshot of the unsteady CFD solution. It shows the relatively high levels of 
eddy viscosity (nondimensionalized by laminar viscosity) in the wake. Mach number contours and vorticity 
magnitude contours (nondimensionalized by speed of sound and unit length) do not indicate any small-scale 
wake structures, wake unsteadiness, or significant wake vorticity. Thus, it is hypothesized here that the 
lower overall lift enhancement predicted by URANS is a direct result of its inability to resolve the wake 
vorticity. Scale-resolving CFD, such as three-dimensional hybrid RANS/LES methods, in combination with 
significantly finer grid resolution in the wake, would be required to improve the results. 

C. Comparison of the unsteady lift of a NACA 0018 in oscillating pitching motion 

The second test case considered here is an airfoil pitching at small angles of attack at a constant free stream. 
The sinusoidal angle of attack profile is a(<f>) = 2° + 2° sin(</>) at low k = 0.08 and medium k = 0.263, 
with Re = 3 • 10 5 . The corresponding theoretical unsteady lift is given in equation (2). Figure 10 shows a 
comparison of a URANS simulation, experiments and the theory of Theodorsen. 2 In the case of low k = 0.08, 
the slope of the CFD curve is still close to 2-7T and thus, close to the quasi steady baseline. This behavior 
suggests that the quasi steady variations are too dominant compared to the unsteady lift effects. The theory 
predicts a reduction of the lift curve slope, which is confirmed by the experimental data. The hysteresis loop 
starts at a = 2° in positive direction towards a = 4°. The loops end at approximately a = 1.9°, visualized 
by the small gap. Although the hysteresis amplitude is similar in all three approaches, the CFD deviates 
significantly. For example, the CFD reveals a clockwise rotational sense. On the contrary, the experiments 
and the theory predict a counter-clockwise rotational sense, which is in agreement with the literature. 4,68 
For the higher k = 0.263, all three approaches yield the same clockwise rotational sense. In theory, the 
switch occurs at approximately k = 0.16. 68 All three curves are tilted more horizontally compared to the 
low k case. Nevertheless, the CFD still shows the highest slope, which is still close to 27r, suggesting that 
the CFD behavior is still dominated by the quasi steady changes. On the other hand, the trend of increased 
hysteresis with increasing k is reproduced by the CFD. Once again, the deficiencies in the URANS may be 
due to its smearing of the wake vorticity and hence not resolving the unsteady wake structures. 

D. Combined free stream oscillation and pitching motion of a NACA 0018 

In the two previous sections, a pure free stream oscillation and a pure pitching motion were investigated. 
Now, pitching motion and free stream oscillations occur in-phase, simultaneously and synchronized. Here, 
both experimental and computational approaches possess some uncertainties, thus, the theoretical unsteady 
lift is considered as well, using van der Wall’s equation (6) based on potential flow theory. This equation 
delineates a non-linear interaction of the two separated effects. In fact, a simple superposition does not 
capture the flow physics correctly. To show this, the NACA 0018 is exposed to a sinusoidal free stream 
«(</>) = tt s (l + 0.5sin(((>)) and an in-phase pitching motion a((f>) = 2° + 2°sin(0), both individually and 
collectively. The time averaged free stream velocity is u s = 13.3 m/s which corresponds to Re = 3 • 10 5 . The 
reduced frequency is k = 0.099. 

Results of the unsteady lift per unit span are shown in figure 11, which compares experiment, theory, 
and CFD. This plot shows the pure free stream oscillation results at fixed a = 2°, the pure pitching motion 
results at 0° < a < 4°, and the combination. During the pure free stream oscillation (green lines and 
symbols), the CFD underestimates the lift between about 0° < </> < 130° and again beyond (j> > 300°. For 
130° < 4> < 300°, the CFD results are close to the experiment and the theory. Regarding the pure pitching 
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Figure 9. CFD eddy viscosity, Mach number, and vorticity magnitude contours in the wake of a NACA 0018 
airfoil at a = 2° subjected to an unsteady free stream oscillation of a = 0.5 and k = 0.074, shown here near 

4 > = ioo°. 
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Figure 10. Experiments, CFD, and theory for a NACA 0018 oscillating between o(</>) = 2° +2° sin(</>) at Re = 310 5 
k = 0.08 (left) and k = 0.263 (right). 
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motion (blue lines and symbols), the experimental results are almost identical to the theoretical prediction 
of Theodorsen. 2 However, similar to the pitching results shown earlier for different fc, CFD overpredicts the 
unsteady effects near maximum lift and underpredicts them near minimum lift (the CFD minimum is close 
to zero). 

Considering the synchronized, combined case (red lines and symbols), the overall shape appears to be 
reasonably well captured by CFD, but again the lift is too close to zero during the second half. Near 
minimum lift, the CFD is closer to the quasi steady theory using simple superposition (gray line). Near 
maximum lift, the CFD is lower than experiment while the van der Wall theory is higher than experiment by 
similar amounts. Both CFD and theory predict significantly lower maximum lift levels than the quasi steady 
simple superposition, in reasonable agreement with experiment. Thus, it appears that the 2-D URANS — 
although able to qualitatively predict the combined effects of surging and pitching — is unable to capture 
the unsteadiness as accurately as the van der Wall theory, particularly in the range of the lowest angles of 
attack. A lift of almost zero at <f> = 270° suggests that the CFD simulation is dominated by the quasi steady 
effects and the induced lift due to the shed wake vortex sheet has only minor influence. On the other hand, 
the theory of van der Wall agrees very well with experiment throughout the lift cycle of the combined case. 



Figure 11. Unsteady lift of a NACA 0018 due to individual and combined/synchronized pitching motion 

q(</>) = 2° + 2° sin (<f>) and oscillating free stream = u s (l + O.5sin(0)), with Re = 3 ■ 10 5 and k = 0.099. 
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E. Dynamic stall of NACA 0018 


In preliminary 2-D CFD studies, attempts were made on unsteady dynamic stall pitching cases using a(<j)) = 
14.5° + 3° sin(0) at Reynolds number Re = 3 • 10 5 . As expected, because even the steady calculations at this 
high of an angle were so far off from experiment (see figure 4), the unsteady pitching results did not have 
any chance of success. Although not shown here, CFD did show hysteresis, but at much higher lift levels 
than experiment, close to the quasi steady levels. 

Given the earlier steady CFD results, we know that two-dimensional RANS is unable to successfully 
compute the flows at high angles of attack in the dynamic stall regime. No doubt three-dimensional com- 
putations will be required in order to account for the side wall influence. Additionally, RANS/URANS 
turbulence models themselves are generally known to be poor for predicting flows with large separated re- 
gions, and they also smear unsteady wake structures that appear to be important in this context. The 
shed wake vortex sheet and its upstream induced velocities is the main driver of the unsteady lift response. 
Considering the three-dimensional CFD results above, the flow on the NACA 0018 wing at a = 14.5° is 
separated over as much as 40% of the upper surface. In this situation, even a three-dimensional URANS 
computation would likely be unsuccessful. Instead, a three-dimensional hybrid RANS-LES computation 
would be required, for example. This type of model, which resolves the larger turbulent scales of motion, 
is known to better capture the physics of massively separated flows, including wake dynamics. However, 
hybrid RANS-LES can be considerably more time consuming than URANS, so it has not been undertaken 
to date for this project. 


VI. Conclusion 

This three-way experimental, computational and theoretical investigation illustrated a number of impor- 
tant aspects of a thick airfoil subjected to pitching and surging. Pre-stall experiments were conducted at 
large surge amplitudes, typically encountered on rotorcraft blades; and this proved to be an excellent basis 
for the validation of Isaacs’ theory. The comparison was satisfactory, but near the maximum lift overshoot, 
the experiments exhibited a high frequency oscillation. This was traced to the formation and shedding of a 
recirculation bubble upstream of the trailing edge that affected the lift by modifying the Kutta condition. 
It was not clear whether this is caused by low Reynolds number effects, the relatively thick airfoil, or their 
combination; future experiments will be aimed at rectifying these two factors. Two-dimensional computa- 
tions underpredicted the overall lift enhancing effect due to surging because the URANS equations tend to 
“smear” the wake; on the contrary, the theory explicitly models the wake vorticity. 

Excellent correspondence was found between experiments and theory for airfoil pitching as well as com- 
bined pitching and surging. In fact, the comparison presented here appears to be the first clear validation of 
van der Wall’s theoretical results. As in the pure surging case, the URANS produced inferior results, pre- 
sumably again because of its inability to resolve unsteady vortical structures in the wake. Nevertheless, the 
computations illustrated the importance of conducting full three-dimensional computations when comparing 
to wind tunnel experiments at high angles of attack, where in particular the tunnel side wall influence needs 
to be included. 
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